Solute Transport across the Lymphatic Vasculature in a Soft Skin Tissue

Simple Summary The lymphatic system plays a crucial role in maintaining fluid and solute balance in biological tissue, which makes the subcutaneous injection a common approach for delivering therapeutic agents through the lymphatic pathway. The transport of drug solutes is regulated by the interstitial fluid pressure. In our study, a time-efficient poroelastic model is developed to mimic pressure relaxation and build-up in the soft skin tissue. This method can accurately and time-efficiently capture the evolution of pressure, which is validated against both the analytical solution and numerical solution of the previous poroelastic model. The increasing porosity and permeability due to deformation alleviate the high pressure caused by the injection. Furthermore, an improved solute transport model is developed to better address the microscopic properties of the lymphatic vessel membrane. The effects of a varying Stokes radius of drug solute and vessel network structures are investigated on lymphatic uptake. Our comprehensive computation model provides a time-efficient and reliable research tool for studying solute transport into the lymphatic system, which can be utilized to support decision-making regarding lymphatic disturbed diseases. Abstract Convective transport of drug solutes in biological tissues is regulated by the interstitial fluid pressure, which plays a crucial role in drug absorption into the lymphatic system through the subcutaneous (SC) injection. In this paper, an approximate continuum poroelasticity model is developed to simulate the pressure evolution in the soft porous tissue during an SC injection. This poroelastic model mimics the deformation of the tissue by introducing the time variation of the interstitial fluid pressure. The advantage of this method lies in its computational time efficiency and simplicity, and it can accurately model the relaxation of pressure. The interstitial fluid pressure obtained using the proposed model is validated against both the analytical and the numerical solution of the poroelastic tissue model. The decreasing elasticity elongates the relaxation time of pressure, and the sensitivity of pressure relaxation to elasticity decreases with the hydraulic permeability, while the increasing porosity and permeability due to deformation alleviate the high pressure. An improved Kedem–Katchalsky model is developed to study solute transport across the lymphatic vessel network, including convection and diffusion in the multi-layered poroelastic tissue with a hybrid discrete-continuum vessel network embedded inside. At last, the effect of different structures of the lymphatic vessel network, such as fractal trees and Voronoi structure, on the lymphatic uptake is investigated. In this paper, we provide a novel and time-efficient computational model for solute transport across the lymphatic vasculature connecting the microscopic properties of the lymphatic vessel membrane to the macroscopic drug absorption.


Introduction
The lymphatic vascular system returns the interstitial fluid to the systemic circulation in the human body, which plays an essential role in homeostasis and immunization in transport [23]. Another approach is simplifying the stress-strain relation to make the computations more efficient. The linear poroelasticity based on small strain theory is often used to simplify the constitutive Equation [24]. In [25], the deformation and strain are obtained by simplifying the poroelastic constitutive equation to simulate the deformation of a hydrogel. In this paper, we develop a new approximate poroelastic model to mimic the relaxation of the pressure in soft tissues.
As one of the two vascular circulatory systems in the human body, the lymphatic vessel network consists of initial lymphatics with only one single layer of lymphatic endothelial cells (LECs), unlike blood vessels which have a complete basement membrane [2][3][4]. For lymphatic vessels, LECs are loosely connected and have a primary valve structure to ensure the one-directional fluid flows from the interstitial space into the lymphatics. In our previous work [13], three different transport conditions are used to investigate solute transport across the lymphatic vessel membrane, including the Kedem-Katchalsky model. In [26,27], an integral form of the modified Kedem-Katchalsky formulation is used to study dermal clearance through the blood and lymphatic vessels. To further address the special physiological structure of the lymphatic vessel membrane, in this paper, we develop an improved Kedem-Katchalsky model for solute transport across the explicit lymphatic vessel wall. The morphology of the initial lymphatic vessel network varies over the location and species [3]. In human skin tissue, blind-ended and interconnected structures (reticular plexus) are observed for initial lymphatic vessels [28][29][30][31]. However, there is no quantitative investigation on the effect of various structures of the lymphatic vessel network on solute transport and absorption through the lymphatic vessel network.
In this paper, an approximate continuum poroelastic model is developed to simulate the pressure evolution in the soft subcutaneous tissue. The pressure evolution of the proposed model is quantitatively validated against the analytical solution and the numerical solution of the poroelasticity model. The effect of elasticity, variable porosity and permeability on pressure build-up is numerically investigated. The convective and diffusive transport of drug solute into the lymphatic vessels is investigated through an improved Kedem-Katchalsky model by implementing the hybrid vessel network in the continuum poroelastic model. At last, heterogeneous vessel networks with various structures are used to investigate the effect of the morphology of the lymphatics on solute transport into the lymphatic system in the multi-layer soft skin tissue.

Methodology
To simulate the microcirculation in the subcutaneous tissue, the physics-based computational model established in this study is combined with physiological models, i.e., transvascular transport models. A continuum approximate poroelasticity model is developed to simulate the interstitial fluid flow in soft subcutaneous tissue. In this section, the governing equations for fluid flow and drug solute transport in the interstitial space and across the vessel network are summarized.

Poroelastic Tissue Model
An approximate poroelasticity model is developed to include the elasticity and compressibility of the skin tissue. This method is both simple and efficient, making it easy to implement. The computational cost is reduced since the constitutive equation to obtain the displacement of a solid skeleton is not solved. At the same time, it can well capture the relaxation of pressure as accurately as the full poroelasticity model [14]. Here, we show the derivation of the approximate poroelasticity model from the Biot poroelasticity theory [32][33][34]. The total stress τ of the poroelastic medium is the sum of the effective stress of the solid skeleton (responsible for the deformation of the solid part of the tissue) and the local fluid pressure as follows: where σ represents the components of effective stress applied to the solid part of the medium, p represents the pore pressure applied to the fluid, and α is the Biot's coefficient.
Based on Biot's poroelasticity theory [32,34], one form of governing equation for the fluid flow in the poroelastic medium is [14,34,35]: where κ is the permeability of the porous medium µ is the viscosity of the fluid, ε v is the volumetric strain and M is the bulk modulus of the fluid-solid mixture.
For the approximate poroelasticity model, we assume that solid displacement is small, the Darcy velocity of the interstitial fluid is small, and the non-zero component of the solid displacement ı is mainly in the normal direction. Based on our homogeneous and small strain assumption (linear elasticity), we have σ = Eε in one dimension, i.e., Hooke's law for the elastic response, where E is Young's modulus of the elasticity. For a 3-D simulation, the volumetric dilation is ε v = ∑ ε i = ∂ı i /∂x i [36], where ı is the displacement of the solid. We have the mean . After applying the stress balance (E v ε v = αp), we obtain the approximation for the time derivative of the volumetric strain ε v as below: Finally, the pressure equation for the poroelastic tissue becomes: The mass balance in the interstitial space with solid displacement ı reads as [21], wherei is the time derivative of deformationi = ∂ı/∂t, i.e., the velocity of the solid skeleton u s . We define 1/D e = α/E v + 1/M, where D e indicates the general softness of the tissue. After combining Equations (5) and (6), we have the governing equation for the approximate deformation: The details of the injection model and the qualification of the proposed continuum poroelasticity model can be found in Appendices A.1 and A.2, respectively. When utilizing 96-processors parallel computation, it takes approximately one day to complete a simulation using the approximate poroelasticity model, which models the transport and pressure evolution in a poroelastic tissue up to 70 s after the start of the injection. The efficiency of our proposed model is comparable to Darcy's law model for the porous medium [11]. However, it takes three days to finish the same simulation solving the constitutive equation of poroelasticity [12]. The approximate poroelasticity model improves the computational efficiency by three times compared to the full poroelasticity model.

Transport Equations
A comprehensive model for the extravascular transport of fluid and drug solutes and the transvascular exchange in the tissue is described for both the continuum medium over the length scale of O (1 cm) and the discrete vessels over the length of O (100 µm) in this Section [11,13].

Fluid Flow in the Interstitial Space
First, Darcy's law is used to describe the fluid flow in the porous media: where u D is the average velocity of interstitial fluid flow in the porous medium, and u D = u f − u s , where u f is the fluid velocity in the pore. On the right-hand side (RHS) of Equation (8), p is the interstitial pressure and µ is the dynamic viscosity of the interstitial fluid. With the flux gain from blood vessels and drainage from lymphatics included, the continuity equation considering the poroelasticity gives [1,[37][38][39][40]: where J v and J l are the net fluid flux gain from blood vessels and net fluid flux loss to lymphatic vessels per tissue unit volume, respectively, and q r is the source term, representing the injection. For the fluid flow in the proposed approximate poroelastic equation, it satisfies Equations (5) and (6), i.e., ∇ ·i = 1/D e ∂p/∂t. The fluid conservation equation for a poroelastic medium without any source term could be written as [5,41]. After substituting Equation (8) and the expression for ∇ ·i into Equation (9), we have the governing equation for pressure: The net fluid fluxes through the blood and lymphatic vessels are determined by Starling's law, which read [37][38][39][40]: where L pb and L pl are the hydraulic conductivity of blood and lymphatic vessels, respectively. Additionally, S/V denotes the surface area of vessels to the unit volume of tissue. Pressures in blood, lymphatics, and the interstitial space are p v , p l , and p i , respectively. In Equation (11), σ is the osmotic reflection coefficient of blood capillary for the solute, where π v and π i are the osmotic pressures of the plasma in the blood vessel and the interstitial space, respectively.

Solute Transport in the Interstitial Space
A general solute transport equation at the macroscopic scale is used for the continuum model [11,13,37,38,42]: where C is the dimensionless drug concentration in the interstitium, normalized by the initial concentration C 0 . The effective molecular diffusion coefficient in the interstitial space is denoted as D. On the right-hand side, φ e denotes the elimination of the drug per tissue unit volume due to binding to the ECM and cells, cell uptake, or chemical reaction, φ b denotes drug perfusion from blood vessels, and φ l denotes drug drainage to lymphatic vessels. Drug absorption through lymphatic vessels is given by φ l = J l C i [38,42,43]. In Equation (13), φ b is set to be zero, considering that large molecules cannot directly transport into blood vessels. The elimination term (φ e ) is set to be zero except in the section studying the role of binding and metabolism to focus on the transport and absorption of free drug molecules.

Transvascular Fluid Flow and Solute Transport
The well-known model for solute transport across the membrane, the Kedem-Katchalsky (K-K) formulation, is used for solute transport across the blood vessels. The fluid and solute fluxes are expressed as: where J B denotes the fluid flux and j bs denotes the solute flux per unit area across the blood capillary wall consisting of the convective and diffusive transport of solute across the vessel wall (i.e., membrane). The convective transport of solute is determined by the convective fluid flow J B | ∂Ω B , while the diffusive transport of solute is determined by the solute permeability P α b . The average concentration across the membrane is calculated by (C + C b )/2, where C b is the concentration of drug solution inside the blood vessels and C is the concentration in the interstitial space. Here, all concentrations are normalized by the initial concentration of the injected drug solution C 0 . The diffusive transport is determined by the difference between the concentration in the interstitial space and the blood vessels ∆C, which is calculated by C − C b . The reflection coefficient σ and the solute permeability P b are used to take into account the permeability of the blood vessel membrane for various solutes. Considering that large molecules (such as monoclonal antibodies with molecular weight over 150 kPa) cannot transport through the blood vessels, we have σ = 1 and P b = 0 to represent a semipermeable membrane [13].
The lymphatic vessel consists of one single layer of endothelial cells lacking the basement membrane and thus, an improved Kedem-Katchalsky (iK-K) model is developed to describe the transport of macromolecules across the lymphatic vessel wall. The fluid flux J L and solute flux j ls across the lymphatic vessel wall are: where δ valve is used to represent the primary valve structure of the initial lymphatics, where slits are opened when the interstitial pressure (ISF) p i is larger than the pressure inside the lymphatics p l . The primary valves openings allow for the transport of large molecules into initial lymphatics, as shown in Figure 1. Solute permeability through the valve P valve is related to the partition coefficient of the solute Φ, which is a function of the Stokes radius of the solute, primary valve opening area fraction δ f , and diffusion across the opening slit D valve .
where δ f is the primary valve opening area fraction, which is the ratio of the area of the valve opening slit and the area of the lymphatic vessel A valve /A vessel per unit length, where Here, w is the opening width, L l is the length of the slit per unit area of the lymphatic vessel wall, and f is the fraction of slit length open to valve width w. Φ(r s ) is the solute partition coefficient, which is a function of the solute Stokes radius, where δ is used to ensure the partition coefficient is positive. δ = 0 when half the width of the opening valve is smaller than the Stokes radius of the solute (w/2 < R s ) and δ = 1 when w/2 > R s . The thickness of the lymphatic membrane is denoted as d l . We use an exponential function to represent the variation of opening width w with the pressure difference between the interstitial pressure and the lymphatics ∆p = p i − p l [26,44].
Where the pressure difference ∆p is in unit kPa. The opening width of the primary valves on the lymphatic vessels w is 10∼40 nm at the basal pressure and the width of the opening slit can reach as high as 64 nm at high interstitial pressure due to the injection according to the simulations in [44]. The baseline values for the parameters in the improved Kedem-Katchalsky model are shown in Table 1 [26,27,44,45].  The tissue and solute properties (such as hydraulic permeability and solute size) and injection process parameters (such as injection duration and volume) are inputs needed for the computational model, as shown in Figure 2. The fluid flow in soft porous tissue is modeled using the proposed poroelasticity model. The solute transport is modeled by the hybrid vessel network model. In the hybrid vessel network model, the diffusion-convection-reaction equation is used to simulate solute transport through the continuum vessels. On the other hand, the improved Kedem-Katchalsky model is developed to simulate the solute transport across the lymphatic vessel membrane for the discrete vessels. With these physics-based computational models as shown in Figure 2, we can predict the absorption rate and evolution of the interstitial fluid pressure [48].

Lymphatic uptake
Absorption rate

Pressure evolution
Physics-based computational models Output SC injection

Solute transport
Hybrid vessel network model Figure 2. Conceptual overview of the physics-based computational models for lymphatic uptake in soft skin tissue through subcutaneous injection.

Computational Setup
The computational domain for the continuum model without explicit vessel network is 80 × 80 × 50 mm, while the computational domain in the case with explicit vessel network is 50 × 50 × 50 mm, including epidermis, dermis, subcutaneous tissue, and muscle layers [49]. The 80 × 80 × 50 mm domain is set to be a symmetric one to allow for a large injection volume (2 mL), where only a quarter of the whole injection site is simulated, with a symmetric boundary condition on the symmetry plane. The injection point P 0 (center of the injection bolus) is located 4 mm below the top surface along the symmetry axis in the single-layered domain, while the injection depth is 5 mm in the multi-layered domain. Unless otherwise specified, the injection volume is 1 mL for all simulations in this paper. The baseline values for all physical and physiological parameters used in the computational model are listed in Table 2, which are of the same order of magnitude as the previous studies [11,13,38,42,[50][51][52], and the units are all SI base units.
The mesh is refined near the injection point and the smallest mesh size used is h = 0.1R, where R is the characteristic length, equal to the radius of the injection bolus, and time step ∆t = 0.001 s. Based on the mesh convergence study, we choose the second good mesh (256 × 256 × 256) and after the local refinement near the explicit vessel walls, the total number of cells is around 31 million, through which we can accurately and efficiently resolve the solute transport across the lymphatic vessels.
The governing equations and computational models are implemented and solved based on the open-source Computational Fluid Dynamics platform, OpenFOAM. The boundary conditions and numerical schemes for the numerical simulation can be found in Appendix A. The key inputs for the model are the tissue and solute transport properties as well as injection process parameters: hydraulic permeability of the porous skin tissue k, elasticity modulus of the soft tissue E, the diffusivity of the drug solute in the ECM D and across the vessel membrane D valve , vascular permeability of blood capillaries L pb , vascular permeability of lymphatic capillaries L pl , the Stokes radius of the drug solute r s , implicit vessel network density S/V, injection volume, injection duration, and injection depth. The values of these parameters are listed in Tables 1 and 2.

Pressure Build-Up and Relaxation
Interstitial fluid pressure (IFP) plays a crucial role in regulating fluid flow and mass transport in the SC tissue. The pressure build-up process in the soft subcutaneous tissue is important to understand the mechanism of drug transport and absorption. The fluid flow and solute transport across the lymphatic vessels and blood vessels are governed by the pressure in the interstitial tissue and vessels, as shown in Equations (11) and (12). Here, the deformable subcutaneous tissue is expressed as a poroelastic medium using the approximate poroelasticity model. In this section, we investigate the effect of various elasticity and permeability, variable permeability, and the heterogeneity of elasticity and permeability in the multi-layered tissue on the pressure build-up during and after the injection process. Absorption and hybrid vessels are not included in these cases.

Elasticity and Permeability
First, we study the effect of elasticity on the development of pressure in subcutaneous tissue. The young's modulus of skin is affected by many factors, such as age, gender, hydration, and skin thickness [53]. The young's modulus or elastic modulus E varies a lot for different layers of the tissue, as shown in Table 3. The young's modulus E of subcutaneous tissue, which is also known as the hypodermis, is around 2 kPa according to [54]. Here, we study how the elastic modulus E affects the pressure build-up in soft and deformable tissue compared to rigid porous tissue modeled by Darcy's law over a wide range of E, assuming the isotropic elasticity in the subcutaneous tissue. From Figure 3, we can see that with smaller elasticity E, the relaxation time is larger, which is proportional to 1/E, and it takes a longer time for the pressure to reach zero and has a lower peak. For a rigid porous medium (E → ∞), which is modeled by Darcy's law, the pressure drops immediately at the end of the injection process t = 5 s. The same sudden decrease happens for the velocity. As the material gets softer with smaller elasticity, the material absorbs the kinematic energy leading to lower pressure and velocity and elongated relaxation time. The interstitial fluid pressure relaxes to a small value, as both solid and fluid phases are incompressible, and the dilatation only happens during the injection [21]. For tissues from different parts of the human body and different animals, the hydraulic conductivity K and the hydraulic permeability k (K = k/µ) vary over a wide range [11,55,[57][58][59]. The magnitude of permeability k varies over a wide range from 10 −16 m 2 to 10 −13 m 2 in the literature. Pressure evolution at three points, i.e., P 0 , P 1 and P 2 along the symmetric axis of the domain, denoted as p 0 , p 1 and p 2 , respectively, are compared for different permeability values in the poroelastic tissue, as shown in Figure 4. With smaller permeability, the pressure peak at point P 0 increases dramatically from about 200 kPa (k = 1 × 10 −13 m 2 ) to 2000 kPa (k = 1 × 10 −14 m 2 ), and for k = 1 × 10 −15 m 2 and k = 1 × 10 −16 m 2 , the peak pressures are 10 Mpa and 14.5 Mpa, respectively, as shown in Figure 4a. Although the magnitude of pressure for different permeability significantly changes, the pressure at the injection point reaches the maximum value almost at the same time, i.e., at the end of the injection. With the decrease of permeability k and increase of the magnitude of pressure, the relaxation time increases dramatically at fixed elasticity E = 80 kPa. For k = 1 × 10 −13 m 2 , the pressure p 0 relaxes to zero quickly after the injection, while p 0 for the tissue with k = 1 × 10 −16 m 2 is still large at the end of observation time t = 20 s. For the other two points, the pressure profile is different for different permeability values. As shown in Figure 4b, the pressure increases at a lower rate P 1 for k = 1 × 10 −14 m 2 , and much lower for k = 1 × 10 −15 m 2 compared to k = 1 × 10 −13 m 2 . At point P 2 , the pressure for k = 1 × 10 −15 m 2 remains almost zero within 20 s, as shown in figure Figure 4c. This is because the fluid can not spread as much farther from the injection point as the permeability decreases to a very small value. Due to the accumulation of fluid, the interstitial pressure P 0 is larger at the injection point for low permeability values.
Hydraulic permeability plays a key role in the magnitude of the pressure. Specifically, the tissue deformation is related to the hydraulic permeability of the porous medium as well as the viscosity of the fluid. When the permeability is rather small, the magnitude of pressure is very large, and a large elastic modulus of the same magnitude as the pressure can significantly affect the pressure evolution. For instance, for k = 1 × 10 −16 m 2 , the peak pressure is about the order of O (1 MPa) when the elastic modulus E reaches O (1 MPa). The evolution of pressure shows a prominent difference compared to rigid porous media with E → ∞, and the relaxation time increases dramatically with decreasing elasticity. When the hydraulic permeability is larger, the pressure magnitude is smaller and it is more sensitive to small elasticity values of the tissue. For instance, for the case with k = 1 × 10 −13 m 2 , the pressure profile changes dramatically when the elastic modulus is as small as O (1 kPa) while it shows little difference compared to rigid porous media when the elastic modulus is of order O (1 MPa).

Variable Porosity and Permeability
The deformation of soft tissue leads to a variable fluid volume fraction or porosity in the tissue, which can be described by Equation (21). Porosity ε, here, is not uniform and constant in this case, and it varies depending on the tissue deformation. The evolution of porosity ε could be described by Equation (21), as shown in the following equation: where we have the gradient of deformation ∇D = (1/M + α/E v )p for the proposed poroelastic tissue model. The hydraulic permeability of the tissue depends on the properties of the tissue matrix, including the porosity and fiber spacing [57]. Various forms of permeability k dependence on the strain have been proposed in the literature [60]. Here, we adopt the Kozeny-Carman (equation Equation (22)) to account for the spatially variable permeability as a function of porosity, which intrinsically depends on pressure. Variable porosity and permeability have been assumed to be a function of pressure, which should be experimentally determined [61]. The domain size used here is the same as the validation section, which is 80 × 80 × 100 mm, and we focus on the effect of porosity-dependent permeability.
The time evolution of the pressure p and Darcy velocity u D are compared for porositydependent permeability and constant permeability values, as shown in Figure 5. The initial value of porosity ε 0 is set to be 0.01, and permeability k 0 varies in a range from 1 × 10 −15 m 2 to 1 × 10 −15 m 2 , respectively. For a constant permeability case, k = 1 × 10 −13 m 2 . During the injection process, permeability changes with pressure, following Equation (21), and the hydraulic permeability k of the porous medium varies following Equation (22). As shown in Figure 6a,b, the magnitude of the maximum pressure for variable permeability with k 0 = 1 × 10 −13 m 2 is less than half that of constant permeability. The velocity behaves the opposite; the maximum velocity of variable permeability is twice that of constant porosity and permeability. With smaller initial permeability, the pressure peak increases, and velocity decreases. Especially for k = 1 × 10 −15 m 2 , the pressure peak rapidly reaches a much higher value. The distribution of permeability and porosity in a single-layered domain is displayed to illustrate the effect of the injection. As shown in Figure 6, with increasing pressure due to the injection of fluid, the porosity increase due to the deformation of the tissue. The increasing porosity leads to a higher permeability. Porosity ε and permeability k are much larger near the injection point, as shown in Figure 6b,c. The increasing permeability helps mitigate the large fluid pressure at the injection point, where fluids could permeate through the tissue and spread out more easily.
The interstitial fluid pressure inside the tissues can be easily affected by the measurement and experimental equipment [17,62,63]. There are very limited quantitative measurements or biological data on the dynamic variation of pressure [17]. The evolution of pressure obtained from the in vivo porcine tissue injection in [17] shows similar profiles as in our simulations.

Transport and Lymphatic Uptake of Drug
In this section, the effect of the heterogeneity of tissue properties in different layers and the heterogeneous vessel network with different structures on the transport and absorption of drug molecules are investigated in poroelastic tissue with variable porosity and permeability. Starting from this section, the drug absorption and the hybrid vessel networks are included.

Multi-Layered Tissue with an Implicit Vessel Network
A multi-layered computational domain is used in this section to make the simulation closer to the real tissue. As shown in Table 3, the elasticity for different skin layers varies in a wide range. In our model, the elasticity for the epidermis, dermis, subcutaneous tissue, and muscle layers is set to be 35 kPa, 10 kPa, and 80 kPa, respectively. Based on a previous study on the thickness of skin layers [49], the thickness of different skin layers varies by many factors, such as the location in the body, gender, and hydration. Here, we use the mean values of the abdomen and thigh for each skin layer. The thickness of the dermis and subcutaneous layers is 2.2 mm and 13.9 mm, respectively, and the muscle layer thickness is set to be 19.9 mm. The hydraulic permeability for each layer varies from κ 0 = 1 × 10 −16 m 2 for the dermis layer to κ 0 = 1.0 × 10 −13 m 2 for the subcutaneous tissue layer [57], as listed in Table 4. We first investigate the drug transport and the lymphatic uptake through the implicit vessel network in a multi-layered poroelastic domain, where surface area per unit volume S/V (implicit vessel network) varies in different layers. The layer-specific mechanical properties affect the pressure build-up and relaxation in soft tissues [12,14,20]. To understand the convective transport driven by the pressure in the interstitial space for the continuum model, the spatial distribution of interstitial fluid pressure, Reynolds number, and the absorption due to lymphatic drainage are displayed at different moments in Figure 7. The interstitial pressure peaks at the injection point and the magnitude of interstitial pressure decreases gradually with time, as shown in Figure 7a. The Péclet number (Pe) represents the ratio of convective and diffusive transport. The distribution of the Péclet number in the domain is used to indicate the convective transport driven by the interstitial pressure (ISF). The convective transport from the interstitial space into the lymphatic vessels is stronger due to the large interstitial pressure during the injection. At the end of the injection t = 5 s, Pe reaches O (10 4 ), where the convection is much stronger than the diffusion in the continuum domain. After the injection, the convection gradually decreases as the interstitial pressure decreases to O (100 Pa) at t = 20 s and t = 60 s. The convective transport of fluid under pressure difference leads to lymphatic uptake. The lymphatic drainage is strong under high interstitial pressure; however, the injection duration is short. As shown in Figure 7c, the magnitude of lymphatic drainage decreases from 1 × 10 −5 s −1 at t = 5 s to 1 × 10 −7 s −1 at t = 60 s.

Multi-Layered Tissue with a Hybrid Vessel Network
The injection leads to increased interstitial pressure which gradually relaxes after the injection, as shown in Figures 6 and 7 in previous sections. The pressure difference between the interstitial space and the lymphatic vessel leads to a convective fluid flow into the lymphatic vessels. The evolution of interstitial fluid pressure is affected by mechanical properties, such as elasticity and permeability, as shown in previous sections. In this section, we study the convective transport and absorption of drug molecules through a 3D hybrid discrete-continuum vessel network in multi-layered soft skin tissue. The details about the generation of the explicit vessel network and the corresponding computational setup can be found in our previous papers [11,13].

Improved Kedem-Katchalsky Model
To illustrate the advantage of the improved Kedem-Katchalsky model in addressing the microscopic properties of the lymphatic membrane, we investigate the convection and diffusion of various molecules with different molecular weights and diffusivity. The hydrodynamic radius of a monomeric immunoglobulin with a molecular weight of 15 kDa is normally around 5 to 6 nm [64]. For albumin with a molecular weight of 66 kDa, the Stokes radius is 3.3 nm [47], and the solute diffusivity is around 5 × 10 −11 m 2 /s [65]. For proteins with molecular weight around 1000 kDa, the Stokes radius is around 9 nm. In this section, 0.2 mL drug solution is injected into the cubic domain (50 × 50 × 50 mm). A discrete vessel network for blood vessels and lymphatic vessels with uniform radius r v = 80 nm, as shown in Figure 8, is embedded into the multi-layered domain. We investigate the transport of various drug solutes by varying Stokes radius r s from 3 nm to 9 nm and the solute diffusivity across the lymphatic valve D valve from 5 × 10 −11 to 1 × 10 −12 m 2 /s. The solute diffusion across the lymphatic endothelial cells is significantly affected by the convective fluid flow into the lymphatics because the opening of the primary valves on the lymphatic vessel wall is determined by the pressure difference, i.e., the convective flow. The diffusion through the opening slits on the lymphatic vessels is also affected the convective lymphatic flow inside the lumen (rapid lymphatic drainage), and thus, it is called convective diffusion here.
The convective and diffusive transport across the lymphatic vessel is analyzed here through the local Péclet number, Pe l . Here, Pe l is defined as U l L/D valve = U l /P valve , where U l is the velocity across the lymphatic vessel wall. The magnitude of U l can be evaluated as U l ∼ εL pl ∆p, which is O (10 −10 m/s) at the end of the injection. The solute diffusivity across the lymphatic vessel D valve divided by characteristic length scale L can be evaluated by the solute permeability P valve , the magnitude of which is O (10 −9 m/s) for the baseline case with D valve = 10 −11 m 2 /s. Thus, the local Péclet number is 1 for the baseline case representing the diffusive and convective transport of solute are of the same magnitude. For the solute diffusivity D valve = 10 −12 m 2 /s, the ratio of convective and diffusive transport is 1.
The evolution of normalized drug volume, i.e., drug percentage, in the tissue for various molecules with different Stokes radii r s and solute diffusivity D s is compared in Figure 9a. The absorption rate is larger for drugs with smaller Stokes radii and larger diffusivity. The characteristic exponential decay time τ d (i.e., time reaches characteristic time for exponential decay, drug volume reaches 36.8% of the initial volume V = e −1 V 0 = 36.8% V 0 ) within 300 s is used to indicate the drug absorption rates [13], which are 16 h, 51 h, and 98 h for Rs = 3 nm with D s = 1 × 10 −10 m 2 /s, Rs = 5 nm with D s = 1 × 10 −11 m 2 /s and Rs = 9 nm with D s = 1 × 10 −12 m 2 /s, respectively. At the same convective flow, the pressure difference is the same between the interstitial space and the lymphatics. However, the solute permeability for smaller solutes is higher as the larger partition coefficient due to the opening of the primary valves and the solute diffusivity is larger, leading to a higher absorption rate. To illustrate the effects of the larger partition coefficient and diffusivity due to decreasing Stokes radii separately, the drug characteristic exponential decay time within 10 s for various solute sizes with fixed diffusivity and different diffusivity are compared in Table 5. In the short term, the effect of solute diffusivity coming from the change of the Stokes Radius is more prominent than the pure change of the Stokes Radius (i.e., the change of partition coefficient). For Rs = 9 nm, τ d increases to 10.3 h compared to Rs = 5 nm with the same diffusivity, causing 18% difference in absorption rate. With decreasing diffusivity D s = 1 × 10 −12 m 2 /s for Rs = 9 nm, τ d increases to 48.8 h, causing 81% difference in absorption rate. In the longer term, the effect of solute diffusivity is smaller as the diffusion slows down due to decreasing concentration gradient and interstitial fluid pressure.  Table 5. Comparison of exponential decay time τ d within 10s between solutes with different Stokes radius and diffusivity. The 'different diffusivity D s ' means that for Rs = 3 nm, Rs = 5 nm and Rs = 9 nm, the solute diffusivity is 1 × 10 −10 , 1 × 10 −11 , and 1 × 10 −12 m 2 /s, respectively. Compared to the Kedem-Katchalsky model with only solute permeability P α l determining the diffusion across the lymphatic vessel wall, the improved Kedem-Katchalsky model can take into account the effects from properties of the vessel membrane and solute molecules on the convective diffusion across the opening of the primary valves, such as the solute size and solute diffusivity across the lymphatic opening slits. From Figure 9a, the absorption rate is lower for solute flux governed by the improved Kedem-Katchalsky model than by the Kedem-Katchalsky model with p l = 1 × 10 −8 m/s.

Molecule Size
In the soft interstitial tissue, the pressure gradually decreases after the end of injection (t = 5 s), as shown in Figure 9b. The time evolution of the interstitial pressure at different points in the domain decreases to around 0.12 kPa at t = 20 s. The interstitial pressure relaxes to the resting pressure p ∞ (82.4 Pa) at around 200 s. Thus, the convective transport and the effect of large interstitial pressure due to the injection in soft skin tissue mainly exist during this period.

Structure of the Lymphatic Vessel Network
The morphology of the lymphatics varies from fractal trees to interconnected irregular shapes depending on the location of biological tissues and species, as reported in [2,3,28,30,66]. In this section, a discrete vessel network for the lymphatic vessels with uniform radius r v = 200 nm is embedded into a domain of size 50 × 50 × 36 mm. Blood vessels are implicitly included in the model. First, the fractal tree structure is generated using the space-fill fractal algorithm [67] to represent the branches of the lymphatic vessels, as shown in Figure 10a. In addition to fractal trees, we use the Voronoi structure to mimic the interconnected initial lymphatic vessel network structure, as shown in Figure 10b. The fractal tree structure and Voronoi structure of lymphatics vessels are embedded into the 5 × 5 × 3 cm domain. Specially, the height of the Voronoi structure is 1.6 cm assuming such a lymphatic structure exists in the dermis and subcutaneous layer. The equivalent surface area of lymphatic vessels per unit volume of the tissue for these two structures is 188.9 m −1 and 114.8 m −1 , respectively. First, we quantitatively compare the evolution of drug volume in tissue for these two structures. As shown in Figure 11, the absorption is much faster through the Voronoi structure than the fractal trees. At the end of the injection, there is 95.8% of the total drug left in the tissue with Voronoi-structured lymphatic vessel network, while that value is 98.6% for fractal trees. During the 5 s injection, the difference in absorption percentage is around 2.8% for these two structures. The characteristic time within 300 s is 18 h and 10 h, respectively, for these two structures. The absorption through the Voronoi-structured vessel network is significantly faster than the fractal tree structure. To further illustrate the high absorption efficiency of the Voronoi structure of the lymphatic network, we investigate the drug transport across the explicit vessel wall. The drug distribution is compared in multi-layered tissue with fractal tree-structured and Voronoi-structured lymphatic vessel networks at t = 100 s. As shown in Figure 12, there are more vessels in contact with the drug plume on the central plane of the domain for the Voronoi-structured lymphatic vessel network. The topology of the Voronoi structure leads to accumulated and concentrated vessels in the dermis and subcutaneous layers.

Conclusions
Interstitial fluid pressure plays a crucial role in governing the convective transport and absorption of drug molecules in the tissue. An approximate continuum poroelasticity model is developed to mimic the pressure build-up and relaxation by introducing the time variation of pressure. The effects of elasticity, porosity, and permeability on pressure evolution are investigated after validating the proposed poroelastic model with the analytical solution and previous computational model. The advantage of this model lies in the simplicity of implementation while maintaining good accuracy in capturing the evolution of the interstitial pressure. The small elasticity of the tissue slows down the relaxation of pressure and lowers the magnitude of pressure, while the magnitude of pressure increases dramatically with decreasing hydraulic permeability. Pressure evolution for low hydraulic permeability tends to be more sensitive to small elasticity, i.e., the relaxation time is longer for low hydraulic permeability at the same elasticity. The effect of variable porosity-dependent permeability is found to alleviate the high pressure due to increased porosity and permeability near the injection site.
To better address the microscopic structure of the lymphatic vessel membrane, an improved Kedem-Katchalsky model is developed for solute transport across the explicit vessel wall. The transport and lymphatic uptake of drug solutes in a multi-layered medium are investigated to elucidate the effect of heterogeneity of multi-layered structure and the heterogeneous vessel network. We investigate the convective diffusion of various macromolecules through the hybrid lymphatic vessel network by varying the Stokes radius and solute diffusivity across the opening lymphatic valves. With a smaller Stokes radius and larger solute diffusivity, the lymphatic uptake is significantly fast. Finally, the effect of the lymphatic vessel network structure is investigated through the fractal trees and Voronoi structures. We find that the Voronoi structure in the subcutaneous layer is more efficient in drug absorption because there are more interactions between the drug molecules and the vessel surfaces in the dermis and hypodermis layers.
For future work, the computational model can be further coupled with the mathematical model for the lymph flow inside the lymphatic vessels. In this way, the effect of the lymph flow rate can be incorporated into the model. The computational models developed here can also be coupled with the pharmacokinetic (PK) models to better predict drug bioavailability.  Institutional Review Board Statement: Not applicable as our study does not involve humans or animals.

Informed Consent Statement: Not applicable as our study does not involve humans.
Data Availability Statement: The data that have been used are confidential.

Acknowledgments:
We thank Mario De Lucio for providing the raw data of pressure in a poroelastic medium and Xiaoxu Zhong for helpful advice on the writing of the initial manuscript. As for the discretization scheme, first, the Euler scheme is applied for the time derivative, and the 'Gauss linear' scheme is applied for the gradient of all fields. For the divergence term, the 'bounded Gauss upwind' is applied for the convection term in the momentum equation, while the 'Gauss limited Vanleer' scheme is applied for the convection term in the transport equation. The boundary of the explicit vessel wall on the side of the interstitial space is denoted as ∂Ω B for blood vessels and ∂Ω L for the lymphatic vessels. The boundary of the domain is denoted as ∂Ω. Variable n is the outside unit normal vector of the boundary surface, pointing to the interstitial space. First, we list the boundary conditions (B.C.s) on the lymphatic vessels. The pressure on the lymphatic vessel wall can be calculated by the fluid flux determined by Starling's law, as shown in Equation (A1) [68]:

Conflicts of
The boundary conditions for the concentration on the lymphatic walls are determined by the improved Kedem-Katchalsky model, i.e., J ls | ∂Ω L = J l (1 − σ)(C l + C)/2 + P α l (C − C l ). The fluid flux is determined by Starling's law, i.e., Equation (A1). The boundary condition is shown in Equation (A2) (a), which represents that the diffusion and convection of drug molecules through the vessel wall being equal to the solute flux through the vessel wall, which is governed by the proposed improved K-K model in the methodology section: The boundary conditions on the blood vessel walls are listed below. The fluid velocity on the blood vessel wall is modeled by Starling's law: similar to the boundary condition for pressure on the lymphatic vessel, the pressure on the blood vessel is as follows: To address the fluid balance after the injection in the long term, we assume J b = J l , which gives a non-zero resting interstitial pressure p ∞ . The equivalence of lymphatic drainage rate and blood perfusion rate is shown in Equation (A5), substituting the values for all parameters and solving the above fluid balance equation, we can obtain a non-zero interstitial pressure at equilibrium p ∞ . For our baseline case, p ∞ = 82.4 Pa.

Appendix A.2. Injection Model
Instead of resolving the needle geometry where the drug is injected, a spatial and temporal variable spherical injection profile is adopted, which represents the spherical fluid bolus [11,13,14]. The injection model with spatial variable injection rate is achieved by the hyperbolic tangent function of distance from the injection point, as shown in (see Equation (A6)). For the injection rate, the smoothed injection rate profile based on hyperbolic tangent functions is adopted here Equation (A7) to help alleviate the sharp change of mass and enhance the convergence of the numerical algorithm.
where the volumetric injection rate q 0 , time duration t i , and total injection volume V i need to be specified. The volumetric injection rate q 0 = V i /t i /( 4 3 πR i 3 ), where R i is the radius of the injection 'sphere'. A previous study [69] showed that longer injection time and larger needle radius help to reduce the injection force and ease patients' pain. For the injection model described by Equation (A6), the sphere represents the fluid plume coming out of the needle, and the radius R i could vary with the needle, which ranges from 0.08-3.8 mm in inner diameter. Here, we take R i = 1 mm as the baseline value for the sphere plume.
where t s = 0.5 and t e = 5.0 is the start and end time for the injection. A and B are adjustable constants, here, we take A = 5/9, B = 7, the same as in [14].

Appendix A.3. Qualification of the Poroelasticity Model
Three observation points on the axis of symmetry are chosen to study the evolution of pressure, P 0 , P 1 and P 2 , where P 0 is the injection point, and P 1 and P 2 are 2 mm and 10 mm below the injection point, respectively. To verify the validity and accuracy of the proposed approximate poroelasticity model, the results for pressure evolution are quantitatively compared to the analytical solution and numerical results of the full poroelasticity model [14]. Appendix A. 3

.1. Comparison against the Analytical Solution
An analytical solution of Equation (5) is obtained in an infinite domain, along with a homogeneous initial condition, i.e., p(x, t = 0) = 0, using the Fourier transform on space and Laplace transform on time sequentially. For a clear presentation, the hydraulic conductivity K = k/µ is used and injection source term q r is modeled as q d s H h (t i − t)δ h (x − x i ), where q d s represents the injection volume, and H h and δ h are approximations of the Heaviside function and delta function, respectively. Here, the Fourier transform of a variable f is denoted bỹ , where ξ is the independent variable in the Fourier space, and its Laplace transform is denoted byf (s) = L[ f (t)], where s is the independent variable in the Laplace space. After applying the Fourier transform to Equation (5), and its initial condition, and then the Laplace transform, the following algebraic equation forp is obtained Immediately,p is solved from the above equation, namely: The Fourier transform of p, i.e.,p, is obtained by performing the inverse Laplace transform top. Noticing thatp is a product of the Laplace transforms of two functions, the convolution theorem is applied andp is Similarly,p is also a product of the Fourier transforms of two functions. After applying the convolution theorem again, the solution of Equation (5) is obtained, namely: Since integrals in the analytical solution Equation (A11) can rarely be computed exactly, we have to evaluate them numerically, given H h (t i − t) and δ h (x − x i ). The spatial integral is performed in a spherical coordinate, namely: 4De K(t−τ) r 2 sin(θ)dr, with a given τ, t, x, and x i , and the Gaussian quadrature rule is applied. Here, R I should be infinite theoretically but is finite in practice, because the smoothed delta function δ h (x) is non-zero only in a small region around x. Then, Equation (A11) is rewritten as: where To compute the integral in Equation (A13), the compound mid-point rule is employed, namely, with t given, where ∆t is the time step and N t = t/∆t is the number of partitioned time intervals.
In the present study, we use H h (t i − t) = q t (t) and δ h (x − x i ) = q s (x), which are detailed in Appendix A.2, although the analytical solution Equation (A11) works for any other smoothed or even the exact Heaviside and delta functions. The reported results used 10 Gaussian quadrature points along each coordinate and ∆t was fixed to be 0.05, and we have checked that almost the same results were obtained using half the number of the Gaussian quadrature points and doubled the time step ∆t. Using the same boundary condition as the analytical solution, we obtain the evolution of the pressure for the proposed poroelasticity model at the point P 2 . The pressure evolution shows good agreement with the analytical solution, as shown in Figure A1. To further verify the effectiveness and accuracy of the proposed model, we also quantitatively compare the results of pressure evolution against the previously obtained poroelasticity model. To further validate the accuracy of the approximate poroelasticity model, we quantitatively compare our results to [14], where the full poroelasticity model was implemented. We use the same injection model, a sphere bolus with radius R = 0.1 mm. Here, the domain size and injection volume are 80 × 80 × 100 mm and 2 mL, respectively. Elasticity and permeability are E = 80 kPa, E v = 267 kPa and k = 1 × 10 −13 m 2 , respectively. As shown in Figure A2, the time evolution of pressure at point P 2 , 10 mm below the injection point, using the approximate poroelasticity model, agrees well with the results using the full poroelasticity model. The small discrepancy could also come from the numerical method. The finite element method is used in [14], while in our simulation, the discretization is based on the finite volume method. Figure A2. Numerical simulation results for pressure using approximate poroelasticity compared to results using full poroelasticity using an identical injection model and parameters [14].

Appendix A.4. Non-Darcy Flow Effect
By implementing the Darcy-Brinkman-Forchheimer equation, the non-Darcy flow effect is investigated on the pressure build-up. Darcy's law is simple and widely used to simulate fluid flow across the continuous porous medium. However, it cannot take into account the boundary effect and the inertia effect [70]. The extracellular matrix is full of fibrous proteins, such as collagen and elastin. Significant departures from Darcy's law first occur at Reynolds numbers on the order of O (1) [71] in the fibrous porous medium. The effect of boundary is important for two layers of fluid-porous medium [72] and flow across a microvessel wall [20].
The axial pressure gradient can be represented as the sum of two terms, one is linear in the velocity (viscous contribution), and the second quadratic in the velocity (inertia contribution) [70,73]. The inertia effect is large during the injection, as the pressure and fluid velocity are large. Thus, the inertia drag may not be neglected in such a condition. Here, the Darcy-Brinkman-Forchheimer equation is implemented to investigate the role of viscous and inertial effects on pressure increase and relaxation. Darcy's law is typically used to approximate fluid flow in the porous medium with a small velocity or Reynolds number.
The Brinkman equation is used to replace Darcy's law to take into account the effect of the boundary, the top surface of the domain, and the porous vessel wall. During the injection the velocity of fluids is rather large, i.e., O (1 m/s), and the Darcy velocity reaches O (10 mm/s) at the end of the injection. The Reynolds' number, Re = ρUR µ , can reach to O (10). The advantages and applicable tissue types of the Darcy model, the Darcy-Brinkman model, and the Darcy-Brinkman-Forchheimer model are listed in the table below [70]. Here, we take into account both the inertial and viscous loss. The general equation governing fluid flow, i.e., the Darcy-Brinkman-Forchheimer equation, reads: where, the magnitude of velocity |u D | = √ u D · u D , F is the Forchheimer coefficient, and µ e is the effective viscosity. Balancing the Darcy damping term µ k u D and the Brinkman viscous term µ e ε ∇ 2 u D , we can obtain the length scale L b = √ k/ε. Over this length scale, the viscous effect is important. The darcy number, defined as Da = k/L 2 , is widely used to quantify the ratio of pore scale and the characteristic length scale [72]. Setting the ratio of the effective viscosity and the fluid viscosity as β, β = µ e /µ, β/εDa is the ratio of the Darcy damping term and the Brinkman viscous term. Taking the characteristic length scale L = L b and β = 1, we have β/εDa = 1, which indicates the viscous force balances out the damping force of the porous medium. Taking the thickness of the skin layer O (1 mm) as the characteristic length scale, β/εDa = β10 −5 . For the scaling analysis, the viscous effect over the large length scale is important when the effective viscosity ratio reaches 10 5 . Thus, we vary the viscosity ratio β from 1 to 10 5 . The evolution of pressure and velocity for different viscosity ratios β is shown in Figure A3. The numerical results are consistent with the above theoretical analysis. The effect of viscosity on pressure is negligible until the viscosity ratio β increases to 10 5 . However, the difference is still small for both pressure and velocity at such a large viscosity ratio. Comparing the Darcy damping term and Forchheimer inertial term, we define the Forchheimer number as F r = ρεFk/2µ. When u D is larger than or of the same order of magnitude as 1/F r (i.e., u D ∼ F r |u D |u D ), the Forchheimer term is as important as the Darcy damping, and the inertial effect is not negligible. Varying the Forchheimer coefficient F from 10 9 to 10 13 , the Forchheimer number ranges from 1 to 10 4 . The pressure at point 2 mm below the injection point p 2 varies non-monotonically with the Forchheimer number. Pressure p 2 first increases with F as the inertial drag increases and the fluid accumulates near the injection point. As F increases to 10 13 , the pressure profile changes as the fluid can hardly spread to the observation point. The Forchheimer coefficient F is closely related to the Reynolds number. For the injection in the subcutaneous layer with maximum Reynolds number Re = 1-10, the magnitude of F is 10 4 , which means the inertial effect is negligible. Figure A4. Effect of inertial forces with various Forchheimer coefficient F on the evolution of (a) interstitial fluid pressure p 2 and (b) Magnitude of Darcy velocity |u D | at point P 2 , 10 mm below the injection point.